#include <iostream>
#include <cmath>

const double PI = (std::atan(1.0) * 4.0);

double source(double x, double y);

double u_b(double x, double y);

double source(double x, double y)
{
    return 2.0 * PI * PI * std::sin(PI * x) * std::sin(PI * y);
};

double u_b(double x, double y)
{
    return std::sin(PI * x) * std::sin(PI * y);
};

int main(int argc, char *argv[])
{
    double x0 = 0.0;
    double y0 = 0.0;
    double x1 = 1.0;
    double y1 = 1.0;

    int N = 20;
    double h = (x1 - x0) / (N - 1);
    int m = N - 2;
    int n = m * m;

    double U[m][m];
    double err = 1.0;
    double tol = 1e-12;

    while (err > tol)
    {
	err = 0.0;
	for (int i = 0; i < m; i++)
	    for (int j = 0; j < m; j++)
	    {
		double lx = (j + 1) * h;
		double ly = (i + 1) * h;
		double old_U = U[i][j];	    
		U[i][j] = source(lx, ly) * h * h;
		if (j != 0)
		    U[i][j] += U[i][j - 1];
		else
		    U[i][j] += u_b(lx - h, ly);
		if (j != m - 1)
		    U[i][j] += U[i][j + 1];
		else
		    U[i][j] += u_b(lx + h, ly);
		if (i != 0)
		    U[i][j] += U[i - 1][j];
		else
		    U[i][j] += u_b(lx, ly - h);
		if (i != m - 1)
		    U[i][j] += U[i + 1][j];
		else
		    U[i][j] += u_b(lx, ly + h);
		U[i][j] /= 4.0;
		double e1 = (U[i][j] - old_U);
		err += e1 * e1;
	    }
	err = std::sqrt(err);
    }
    double error = 0.0;
    for (int i = 0; i < m; i++)
    {
	for (int j = 0; j < m; j++)
	{
	    double lx = (j + 1) * h;
	    double ly = (i + 1) * h;
	    double e1 = u_b(lx, ly) - U[i][j];
	    error += e1 * e1;
	}
    }
    error = std::sqrt(error);
    std::cout << "error = " << error << ";" << std::endl;

    std::cout << "U = [" << std::endl;
    for (int i = 0; i < m; i++)
    {
	for (int j = 0; j < m; j++)
	    std::cout << U[i][j] << " ";
	std::cout << std::endl;
    }
    std::cout << "];" << std::endl;;
    return 0;
};
